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Abstract 

Factorization  of  linear  programming  (LP)  models  enables  a  large  portion  of  the  LP 
tableau  to  be  represented  implicitly  and  generated  from  the  remaining  explicit  part.  Dy- 
namic factorization  admits  algebraic  elements  which  change  in  dimension  during  the  course 
of  solution.  A  unifying  mathematical  framework  for  dynamic  row  factorization  is  presented 
with  three  algorithms  which  derive  from  different  LP  model  row  structures:  generalized 
upper  bound  rows,  pure  network  rows,  and  generalized  network  rows.  Each  of  these  struc- 
tures is  a  generalization  of  its  predecessors,  and  each  corresponding  algorithm  exhibits  just 
enough  additional  richness  to  accommodate  the  structure  at  hand  within  the  unified  frame- 
work. Implementation  and  computational  results  are  presented  for  a  variety  of  real-world 
models.  These  results  suggest  that  each  of  these  algorithms  is  superior  to  the  traditional, 
non-factorized  approach,  with  the  degree  of  improvement  depending  upon  the  size  and  qual- 
ity of  the  row  factorization  identified. 


1      Introduction 

A  recurring  theme  in  the  development  of  algorithms  for  linear  programming  has  been  the  iden- 
tification and  exploitation  of  special  problem  structure.  Ideas  as  apparently  disparate  as  the 
bounded- variable  simplex  method,  primal  and  dual  decomposition  methods,  pure  and  gener- 
alized network  primal  simplex  algorithms,  primal  partitioning  and  column  generation  schemes 
may  be  unified  to  a  degree  with  this  view. 

The  factorization  approach  introduced  by  Graves  and  McBride  [1976]  isolates  special  struc- 
ture in  LP  tableaus.  We  are  interested  in  using  factorization  to  reinterpret  existing  algorithms, 
and  to  discover  common  principles  and  apply  them  to  develop  new  algorithms.  Although  all 
algorithms  developed  this  way  will,  in  theory,  solve  any  LP,  the  efficiency  of  any  particular  fac- 
torization approach  will  be  influenced  by  the  relative  number  of  factored  constraints  and  their 
influence  on  the  algorithm:  the  size  and  quality  of  the  special  structure  isolated  determines 
the  influence  of  any  particular  factorization  applied  to  any  particular  LP. 

Based  on  prior  work  by  Brown  and  Graves  [1975],  in  which  generalized  upper  bound  rows 
were  successfully  incorporated  in  a  large-scale  optimization  system,  we  are  interested  in  pursuing 
dynamic  row  factorization,  where  the  dimension  of  the  factored  structure  may  vary  (or  even 
fail  to  be  present)  as  the  solution  progresses.  In  our  setting,  we  require  the  row  structure  of 
the  model  instance  to  be  specified  prior  to  solution,  and  that  this  structure  remain  fixed  during 
solution.  An  extension  of  this  approach  is  to  allow  the  row  structure  to  vary  as  the  model  is 
solved:  this  is  a  conceptually  simple  extension  of  the  approach. 

Each  algorithm  is  developed  by  factoring  the  constraints  of  the  LP  model  into  two  classes: 
those  that  have  the  special  structure  {factored)  and  those  that  do  not  {explicit).  This  constraint 
factorization  induces  a  factored  structure  in  the  LP  tableaus  which  is  exploited  computationally. 
We  demonstrate  the  dynamic  factorization  approach  for  three  special  structures: 

generalized  upper  bound  rows; 

pure  network  rows;  and 

generalized  network  rows. 

We  implement  each  of  the  factorization  algorithms  by  integrating  it  within  the  X-System 
(Brown  and  Graves [1975]). 

While  the  terms  "partitioning"  and  "factorization"  are  frequently  used  interchangeably  in 
the  literature,  we  observe  a  distinction  between  the  two  approaches.  We  consider  partitioning 
methods  to  be  based  on  special  structure  in  the  original  problem  instance,  which  need  not 
induce  special  structure  in  the  LP  tableau — in  fact,  the  method  need  not  be  tableau- based. 
In  contrast,  factorization  methods  are  based  on  special  structure  which  occurs  in  bases  and 
thus  in  the  basic  tableau.  Thus,  we  classify  dual  decomposition  (Dantzig  and  Wolfe  [I960]), 
primal  decomposition  (Benders  [1962]),  and  primal  partitioning  (Rosen  [1964])  as  examples  of 
partitioning  methods. 

Perhaps  the  earliest  example  of  what  we  consider  factorization  is  the  treatment  of  simple 
upper  bounds  by  Dantzig  [1954]  and  [1963]  and,  independently,  by  Charnes  and  Lemke  [1954]. 
They  observe  that  it  is  more  efficient  to  enforce  the  "logical"  upper  bound  constraints  with 
logical  tests  within  the  algorithm  rather  than  treat  them  explicitly  along  with  other  "structural" 
constraints.  While  not  originally  presented  in  the  context  of  a  formal  tableau  factorization,  the 
approach  is  easily  viewed  as  such. 

The  mutual  primal-dual  method  of  Graves  [1965]  focuses  attention  on  the  special  role  of 
nonnegativity  constraints  in  linear  programming.  A  clear  distinction  is  drawn  between  the 
computational  convenience  of  treating  nonnegativity  constraints  implicitly  rather  than  explic- 
itly and  the  unambiguous  mathematical  equivalence  of  all  problem  constraints,  structural  or 
nonnegativity.  Emphasizing  the  special  importance  of  inequality  constraints,  the  approach 
yields  an  elegant  theory  and,  as  we  will  see,  efficient  implementations.  We  view  this  algorithm 
as  the  first  formal  example  of  factorization. 


A  similar  primal-dual  algorithm  is  presented  by  Balinski  and  Gomory  [1965].  Related  work, 
in  which  efforts  are  made  to  exclude  slacks  from  the  product-form  representation  of  the  primal 
basis,  includes  that  of  Zoutendijk  [1970]  and  Powell  [1975]. 

Dantzig  and  Van  Slyke  [1967]  extend  the  earlier  work  for  simple  upper  bounds  and  lend 
a  more  structured  treatment  to  generalized  upper  bounds  (GUB).  In  a  problem  with  p  GUB 
constraints  and  m  structural  constraints,  their  approach  requires  a  working  basis  of  dimension 
(m  -f  1),  a  considerable  savings  when  p  is  large. 

Hartman  and  Lasdon  [1972]  specialize  this  approach  to  the  multicommodity  capacitated 
transshipment  problem.  In  this  case,  the  structure  of  the  basic  pure  network  columns  in- 
troduces additional  structure  into  the  working  basis,  allowing  further  simplifications  in  basis 
representation  and  update  techniques.  Helgason  and  Kennington  [1977]  develop  techniques  for 
representing  the  working  basis  in  product  form  and  provide  graphic  interpretation  of  the  basis 
updates.  Kennington  [1977]  reports  an  implementation  of  the  algorithm. 

Mc Bride  [1972]  and  Graves  and  McBride  [1976]  formalize  and  generalize  the  factorization 
approach.  They  view  factorization  as  a  unifying  framework  for  tableau- based  simplex  special- 
izations and  illustrate  this  by  developing  a  variation  of  the  GUB  algorithm  of  Dantzig  and  Van 
Slyke  and  a  GUB  algorithm  for  the  doubly-coupled  linear  programs  of  Hartman  and  Lasdon 
[1970].  They  present  a  new  algorithm  for  the  set  partitioning  LP  and  an  equality-constrained 
form  of  the  pure  network  with  side  constraints  model.  Brown  and  Graves  [1975]  report  an 
implementation  of  inequality-form,  dynamic  GUB  row  factorization  for  large-scale  problems. 

Schrage  [1975]  extends  the  succession  of  simple  and  generalized  upper  bounds  by  introducing 
variable  upper  bounds  (VUB),  which  are  constraints  of  the  form  x3  <  x^,  where  x^  is  said  to  be 
the  variable  upper  bound  of  xy  Schrage  implicitly  represents  the  VUB  constraints  by  expressing 
VUB  variables  in  terms  of  other  variables.  This  permits  the  basis  representation  to  be  treated  in 
two  parts,  one  a  large  matrix  which  changes  infrequently  and  thus  needs  only  occasional  update, 
and  the  other  a  small  working  basis  which  requires  regular  attention.  Thus,  computation  and 
storage  savings  may  be  realized.  Schrage  [1978]  extends  these  ideas  to  what  he  calls  generalized 
VUB  (GVUB)  constraints,  which  arise  frequently  in  models  with  fixed  charges. 

Klingman  and  Russell  [1975]  sketch  a  factorization  method  for  solving  transportation  prob- 
lems with  side  constraints.  They  suggest  techniques  for  performing  simplex  iterations  and 
updating  the  problem  representation.  Chen  and  Saigal  [1977]  present  a  similar  approach  for 
solving  capacitated  network  flow  problems  with  additional  linear  constraints.  Both  of  these  pre- 
sentations employ  a  graphical  description  of  the  basis  update  and  treat  the  basis  in  two  parts: 
one  corresponding  to  a  rooted  spanning  tree  defined  on  the  underlying  graph,  and  the  other 
a  general  working  basis.  Glover,  et  al.  [1978]  report  an  implementation  of  the  Klingman  and 
Russell  design,  but  one  which  (curiously)  only  accommodates  a  single  side  constraint.  McBride 
[1989]  reports  an  implementation  which  requires  the  pure  network  rows  to  be  equalities  and 
allows  more  than  one  side  constraint. 

Generalized  networks  with  side  constraints  are  addressed  by  Hultz  and  Klingman  [1976], 
who  present  details  for  the  simplex  priceout,  column  generation,  and  basis  update.  Hultz  and 
Klingman  [1978]  report  an  implementation  that  (curiously)  solves  the  "singularly  constrained" 
generalized  network  problem.  McBride  [1989]  reports  an  implementation  that  is  not  restricted 
to  a  single  side  constraint. 

The  factorization  approach  has  been  extended  by  consideration  of  embedded  structures. 
Glover  and  Klingman  [1981]  consider  an  LP  with  embedded  pure  network  structure,  i.e.,  the 
pure  network  structure  appears  in  only  a  subset  of  the  rows  and  columns  of  the  technological 
coefficient  matrix.  They  give  an  algorithm  similar  in  spirit  to  their  pure  network  with  side 
constraint  model,  but  the  presence  of  the  "side  variables"  significantly  complicates  the  basis 
representation  and  update.  They  report  an  implementation  of  the  algorithm  but  (curiously) 
restrict  test  problems  to  have  no  complicating  variables. 

McBride  [1985]  solves  an  LP  with  embedded  generalized  network  structure,  presenting  meth- 
ods for  pricing,  column  generation,  basis  representation  update  and  data  structures.  A  success- 
ful implementation  is  reported  to  be  about  five  times  faster  than  MINOS  (ca.  1977:  Murtagh 
and  Saunders  [1977])  for  the  models  tested. 


Algorithms  to  solve  problems  with  special  substructures  have  motivated  research  to  effi- 
ciently identify  such  substructures.  Brearley,  Mitra  and  Williams  [1975]  describe  algorithms 
for  detecting  GUB  row  sets  and  exclusive- row  structure  sets  (a  set  of  rows  whose  structure 
may  be  transformed  to  GUB  by  column  scaling).  Greenberg  and  Rarick  [1974]  and  Brown  and 
Thomen  [1980]  develop  algorithms  to  identify  GUB  sets.  Brown  and  Wright  [1984]  identify 
pure  network  constraint  substructures.  Brown,  McBride  and  Wood  [1985]  present  a  method  for 
locating  embedded  and  row-only  generalized  network  structures. 

Todd  [1983]  develops  a  geometric  interpretation  of  factorization  which  is  for  our  purposes 
equivalent  to  the  algebraic  development  of  Graves  and  McBride  [1976]. 

In  the  following  sections  we  establish  notational  conventions,  develop  the  mathematical 
foundation  of  a  primal-dual  simplex  method  and  show  the  effects  of  row- factorization.  Next, 
a  general  row- factorization  algorithm  is  developed,  specialized  to  GUB,  pure  network,  and 
generalized  network  rows,  and  tested  on  a  suite  of  real-world  problems.  Finally,  we  discuss  how 
the  methods  can  be  generalized  further  and  how  they  can  be  applied  more  effectively. 

2     Mathematical  Preliminaries 

The  traditional  statement  of  the  linear  programming  (LP)  problem  is 


run  : 

y 

wy 

s.t. 

a,y  <  n 

,  i  =  1,.. 

.  ,m 

e,y>0 

,;  =  !,. 

.  .,n 

(LP) 


where  y  is  an  n- vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  each  a*  an  n-vector 
of  technological  transformation  coefficients,  each  r;  a  scalar  right-hand  side  coefficient,  and  e; 
the  jth  unit  vector.  While  this  statement  of  the  problem  is  clear  and  unambiguous,  there  are 
reasons  for  preferring  an  alternative.  The  insistence  upon  drawing  a  formal  distinction  between 
the  "structural"  constraints  any  <  r%  and  the  "nonnegativity"  constraints  e;y  >  0  obscures 
the  mathematical  structure  of  the  problem  by  suggesting  that  the  two  types  of  constraints  are 
inherently  different.  Certainly  the  exploitation  of  the  special  structure  of  the  e^y  >  0  constraints 
leads  to  computational  efficiencies;  however,  in  our  theoretical  development  of  the  algorithm, 
we  prefer  to  treat  them  simply  as  general  inequality  constraints. 

In  order  to  achieve  a  consistent  form,  we  rewrite  the  nonnegativity  constraints  as  —  e7-y  <  0 
and  group  them  with  the  structural  constraints.  The  problem  statement  then  becomes 


(PLP)     min  :     wy 

s.t.  :     diy  <  Ti 


i  =  1,.  ..,m  +  n, 


where  wy  is  called  the  extremal  function. 

From  the  standpoint  of  a  primal  algorithm,  a  matrix  partitioned  form  of  the  primal  tableau 
is  derived.  Let  {al},  at2, . . . ,  a,n}  be  a  basis  for  Rn  at  y  .  For  notational  convenience  we  will 
partition  the  constraints  into  two  sets,  those  that  are  basic  (binding)  at  y°  and  those  that  are 
nonbasic  (not  necessarily  binding)  at  y° 
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Using  this  notation,  the  current  basic  solution  y°  may  be  expressed  as  By0  —  f,  and  since 
the  rows  of  B  are  by  definition  linearly  independent,  y°  exists. 

To  isolate  the  important  algebraic  components  let  us  assume  that  at  the  current  basic 
solution  the  basis  consists  of  h  structural  constraints  and  (n  —  h)  nonnegativity  constraints. 
Then 
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In  partitioned  matrix  form, 
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and  thus 
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Similarly,  D  can  be  written  as 
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and  in  partitioned  matrix  form 
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and  we  shall  call  DP  the  principal  part  of  the  tableau.  By  partitioning  w  =  (w\,W2),    g    = 
(<?i>  <?2)T  and  y°  =  (y®,  y®),  the  complete  tableau  may  be  written  in  partitioned  matrix  form  as 
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Note  that  yj  is  displayed  explicitly  in  this  tableau.     Also,  y\   —  0  since  the  corresponding 
nonnegativity  constraints  are  basic  and  thus  binding. 
The  corresponding  dual  problem  is 


(DLP)     max  :     xr 

X 

s.t.       xaJ  <  w3 
xei  <  0 


,  j  =  1 , . . .  ,n 
,    i  =  1, . . .  ,m 


To  develop  a  matrix  partitioned  form  of  the  dual  tableau,  we  proceed  as  before.  Assuming 
the  dual  basis  consists  of  h  structural  constraints  and  m  —  h  nonnegativity  constraints,  we  have 

T      =(t\t2,...,tm)      =  (ajl,aj*,...,ajh,eik+\...,eim) 
u     =(ul,v?,...,um)     =(wj\wh,...,wjh,0 0), 

so 


u-  (u\u2)  =  (w\0). 
The  nonbasic  constraints  are  then 

K  =  (k] ,  A:2, . . . ,  kn)  =  (e11 ,  ex\  . . . ,  e\  ajh»  ,...,aj») 


v  =  (v\v%...,vn)  =  (0,0,  ...,0,iu'ft+\...,  «>>»), 


SO  D=  (v1,  V2)  =  (0,  w2)  . 

The  matrix  partitioned  form  of  the  basis  is  then 
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I  \ 

Au        0 

\A2]        I    ) 


}h  , 

}  m  —  h 


Q 


/ 


A; 


\-A2lAul 
with  the  remaining  constraints  forming 


m  —  h 


\ 

0 

1  ) 


}h  , 
}m  —  h  , 


K  = 


The  principal  part  of  the  dual  tableau  is 


/    An 
0     A22 


QK 


*n          0  1 

'  I    AX2' 

-A2]AU]     I  _ 

o    yl22 

Aul                 An1Ai2 

^2Mii       ^22  - 

A 

2\AU  A\2 

which  we  find  to  be  exactly  the  principal  part  of  the  primal  tableau,  so 

DP  =  QK  . 
Thus  a  single  tableau  representation  supports  both  primal  and  dual  algorithms. 

3     Interpretation  of  Primal  and  Dual  Forms 

We  may  interpret  a  primal  or  dual  algorithm  as  simply  different  perspectives  of  this  same 
tableau,  wherein  a  primal  algorithm  basis  change  is  viewed  as  exchanging  primal  constraints  and 
a  dual  basis  change  exchanges  dual  constraints.  The  classical  Simplex  Method  may  then 
be  interpreted  as  solving  (PLP)  using  the  dual  perspective.  That  the  classical  Simplex 
Method  is  naturally  interpreted  as  a  dual  algorithm  comes  as  a  surprise  to  the  conventionally 
trained.  However,  the  consequent  mathematical  insight  is  compelling,  especially  in  light  of 
the  notational  simplification  and  apparent  underlying  role  of  A^ ,  which  we  refer  to  as  the 
transformation  kernel. 

There  are  several  reasons  for  preferring  a  Primal-Dual  algorithm  to  the  Simplex  Method. 
From  a  computational  standpoint,  because  slack  variables  are  carried  logically  rather  than 


introduced  explicitly,  we  are  able  to  clearly  identify  the  essential  information  needed  to  execute 
the  algorithm.  The  matrix  A\l  plays  a  key  role  in  the  calculation  of  the  tableau,  and  the  entire 
tableau  can  be  constructed  from  A\^  and  original  problem  data.  Since  A^^  is  a  submatrix 
of  the  inverse,  T_1,  used  by  the  Simplex  Method,  it  is  smaller  and  requires  fewer  arithmetic 
operations  to  update  than  does  T_1. 

A  second  advantage  of  a  Primal-Dual  Algorithm  lies  in  the  flexibility  it  offers  for  special- 
ization to  particular  problem  classes  or  structures.  Indeed,  it  is  the  special  structure  and 
simplicity  of  the  nonnegativity  constraints  that  motivate  the  development  of  the  algorithm  in 
the  first  place.  It  is  frequently  the  case  that  other  special  structures  can  be  identified  in  classes 
of  (PLP).  Examples  of  such  structures  include  simple  upper  bounds,  generalized  upper  bounds, 
variable  upper  bounds,  pure  and  generalized  network  substructures,  etc.  Such  structure  may  be 
"static"  in  that  its  nature  and  dimension  remains  fixed  throughout  the  solution  process,  or  the 
structure  may  be  "dynamic"  in  which  case  its  precise  nature  and/or  dimension  may  vary  as  the 
problem  is  solved.  Some  special  structures  may  be  more  strongly  characterized  by  their  column 
structure  and  others  by  their  row  structure.  The  Primal-Dual  perspective  leads  naturally  to 
explanations  of  the  implications  of  virtually  any  such  problem  structure  and  greatly  simplifies 
the  implementation  of  such  a  specialization. 

When  an  LP  appears  as  a  subproblem  in  a  more  sophisticated  solution  setting  (for  example, 
in  a  mixed  integer  programming  problem  or  a  nonlinear  programming  problem),  the  row/column 
symmetry  of  a  Primal-Dual  Algorithm  is  of  critical  importance  in  specializing  the  solution 
approach.  The  inherent  symmetry  of  such  an  algorithm  permits  easy  adaptation  to  branch- 
and- bound  and  cutting- plane  approaches  to  mixed  integer  programming,  to  column  generation 
settings,  as  well  as  to  primal  and  dual  decomposition  techniques. 

We  believe  the  reason  for  this  flexibility  offered  by  the  algorithm  lies  in  its  more  complete 
mathematical  foundation.  There  is  a  natural  consistency  that  arises  from  the  choice  of 
a  vector  space  having  the  same  dimension  as  the  problem  variables  that  is  lacking 
in  other  approaches.  A  natural  geometric  interpretation  of  the  solution  trajectory  follows  di- 
rectly from  this  development.  Incidental  issues  such  as  finding  an  initial  basic  feasible 
solution  and  dealing  with  degeneracy  are  resolved  constructively  in  this  math- 
ematical framework  [Graves  1965].  Other  approaches  resort  to  unnecessarily  complicated 
tangential  efforts. 

All  the  research  results  reported  here  can  be  developed,  with  some  effort,  in  the  framework 
of  the  classical  Simplex  Method.  However,  we  choose  to  present  these  results  in  the  manner  of 
their  development  —  the  mutual  Primal-Dual  view  presented  by  Graves. 

4      Column  and  Row  Generation 

Rather  than  maintain  a  complete  tableau  DP,  now  consider  the  generation  of  just  column  c 
of  this  tableau.  Rewriting  in  a  manner  that  highlights  our  intentions,  and  labelling  row  and 
column  partitions  for  identification 


DP 


_  w 


U)  Uj) 


(ii)  \-A2i[Aul]     A22  -  A2i[Aul  An]  J 


By  properly  sequencing  our  computations  we  will  exploit  the  fact  that  region  (ii)  of  a  given 
column  is  simply  a  linear  combination  of  terms  in  region  (i)  of  the  same  column. 

Assume  we  want  to  place  the  current  representation  of  column  c  into  a  work  array  z,  which 
we  partition  as  zT  =  (z{ ,  zj)  to  correspond  to  regions  (i)  and  (ii).  Expressed  in  terms  of  the 
transformation  kernel  A^ ,  we  compute  column  c  as 

if  c  is  in  (j), 


and  then 


or,  if  c  is  in  (jj), 


and  then 


-2  =  -A21  [z\]    ; 
z,  =  [Au\A12y]  , 
z2  =  (A22)c  -  A2\[zi] 


Then  the  current  representation  of  column  c  is  available  in  zT  —  (zf,z%). 

The  computation  of  row  r  of  the  tableau  proceeds  in  a  similar  manner.   We  now  view  the 
principal  part  of  the  tableau  as 


0) 


Uj) 


DP  = 


(«) 


/ 


\ 


\[-^21^n]       ^22+[-^21^nMl2^ 


If  we  want  to  place  the  current  representation  of  row  r  in  a  work  array  5  partitioned 
conformably  with  (j)  and  (jj)  as  £  =  (53,24),  we  compute 
if  row  r  is  in  (i), 


h  =  [^i/lr   , 


and  then 


or,  if  row  r  is  in  («), 


and 


=  \zx\A 
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:3  =  [(-^2iMn1] 


^4   =  (^22  )r  +  foMu     . 


and  the  current  representation  of  row  r  is  available  in  5  =  (z.3,54). 

We  see  that  in  each  case  calculations  proceed  by  first  using  a  representation  of  A^  to 
compute  a  portion  of  the  row  or  column  and  then  using  this  initial  computation  and  original 
problem  data  to  compute  the  remaining  part.  We  will  discover  that  our  specializations  ex- 
tend this  approach  by  introducing  additional  tableau  partitions  which  allow  this  computational 
strategy  to  be  applied  on  a  larger  scale. 


5  Transformation  Kernel  Update 

The  dynamic  behavior  of  A^  is  important.  We  see  from  the  primal  row  basis  B  and  nonbasic 
rows  D  that  the  dimension  of  ^l^1  corresponds  to  the  number  of  basic  structural  constraints, 
or,  equivalently,  to  the  number  of  nonbasic  nonnegativity  constraints  (recall  that  if  a  nonnega- 
tivity  constraint  is  nonbasic  and  thus  nonbinding,  the  corresponding  variable  may  possibly  be 
nonzero).  Recalling  that  our  primal  view  of  a  basis  exchange  is  as  an  exchange  of  constraints 
between  B  and  D,  we  see  that  one  of  four  cases  may  occur  during  a  pivot 

A  structural  constraint  enters  the  basis  B  and  a  structural  constraint  leaves  the  basis  and 
enters  D.  Since  the  number  of  basic  structural  constraints  (and  the  number  of  nonbasic 
nonnegativity  constraints)  remains  unchanged,  the  dimension  of  A^  is  unchanged.  A 
pivot  of  this  type  involves  a  row  in  region  (j)  of  B  and  a  row  in  region  (ii)  of  D,  and  thus 
it  corresponds  to  a  pivot  coordinate  in  the  location  ((ii),  (j))  of  the  tableau  DP. 

A  nonnegativity  constraint  enters  the  basis  and  a  nonnegativity  constraint  leaves  the 
basis.  Again,  the  dimension  of  A^l  remains  unchanged.  Since  this  pivot  involves  a  row 
in  region  (jj)  of  B  and  a  row  in  region  (i)  of  D,  the  corresponding  tableau  DP  pivot 
coordinate  lies  in  ((i),  (ji))- 

A  structural  constraint  enters  the  basis  and  a  nonnegativity  constraint  leaves  the  basis, 
and  thus  the  number  of  basic  structural  constraints  (equivalently,  the  number  of  nonbasic 
nonnegativity  constraints)  increases  by  one.  The  dimension  of  A^  is  increased  by  one. 
This  corresponds  to  a  pivot  coordinate  in  region  ((H),  (jj))  of  the  tableau  DP. 

A  nonnegativity  constraint  enters  the  basis  and  a  structural  constraint  leaves  the  basis, 
and  thus  the  dimension  of  A^1  decreases  by  one.  The  corresponding  pivot  coordinate  in 
DP  is  ((i),  (;)). 

We  see  that  we  may  exert  some  influence  on  the  behavior  of  the  dimension  of  A^]  by  our 
strategy  for  selecting  target  exchanges  for  primal  and  dual  constraints  (i.e.,  our  pricing  strategy) 
and  through  our  tie-breaking  rules  for  choosing  pivot  row/column,  and  that  this  dynamism  is 
an  inherent  feature  of  an  effective  algorithm.  We  have  already  seen  the  fundamental  importance 
of  the  kernel  (An)  in  our  computations.  Thus,  a  successful  implementation  must  manage  this 
dynamic  behavior  efficiently  and  reliably. 

6  Factorization 

The  row-jactorized  problem  to  be  considered  is 

(FLP)     min  :        wy 

s.t.  :        Ey  <  r         }  explicit  constraints 
Fy  <  b        }  factored  constraints 
—  Iy  <  0         }  nonnegativity  constraints      , 

where  y  is  an  n-vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  F  a  matrix  of  con- 
straint coefficients  for  "explicit"  constraints  with  right-hand  side  m-vector  r,  F  a  matrix  of 
constraint  coefficients  for  "factored"  constraints  with  right-hand  side  p- vector  6,  and  -/  the 
negative  of  the  identity  matrix.  In  this  general  development,  we  refer  to  the  F-type  constraints 
as  "factored"  only  to  distinguish  them  from  the  "explicit"  F-type  constraints,  and  assume  noth- 
ing about  their  structure.  Not  until  our  specializations  later  will  we  impose  special  structure  on 
F,  and  the  structures  we  will  consider  may  permit  the  representation  of  the  F-type  constraints 
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without  the  inversion  of  a  matrix.  We  will  see  that  this  approach  is  centered  around  handling 
the  part  of  the  basis  corresponding  to  the  F-type  constraints  explicitly  while  factoring  the  por- 
tion of  the  basis  corresponding  to  the  F-type  constraints.  The  notation  is  chosen  to  suggest 
this  idea. 

Recall  that  a  basis  for  the  primal  algorithm  consists  of  n  linearly  independent  rows  from 
the  constraint  matrix  when  it  is  assumed  to  include  both  structural  (explicit  and  factored)  and 
nonnegativity  constraints.  Assume  that  the  current  row  basis  consists  of  k  rows  from  F,  I  rows 
from  F  and  (n  —  (A:  -I-  /))  rows  from  — /.  Repeating  our  notation 


k+l         n-(k+l) 


B  = 


A} 


\ 

An 

-I     J 


}k  +  l 
}n-(k  +  l) 


where  [An    An]  includes  all  basic  structural  rows,  from  both  E  and  F. 

We  will  ultimately  be  interested  in  isolating  the  effect  of  each  type  of  structural  constraint 
algebraically  in  the  factored  tableau,  and  thus  we  require  greater  resolution  in  our  factored 
basis.  Introducing  obvious  notation,  we  have 


/ 


[Au   An]  = 


l  n-(k+l) 


E\\      En        E\3 

V  F21    F22     E23 

where  the  kernel  of  dimension  (k  +  I)  is  given  by 


}k 


[An]  = 


En 
F2, 


E\2 

F>2 


Because  A\\  is  a  basis  for  Rk+l  it  follows  that,  it  is  always  possible  to  identify  among  the 
columns  of  [F2\  F22]  a  nonsingular  submatrix  F22  of  dimension  /,  since  otherwise  the  rank  of 
[F2\  F22]  is  at  most  (/  -  1)  and  thus  the  rank  of  An  is  at  most  (k  +  I  -  1),  or  equivalently 
the  rank  of  B  is  at  most  (n  -  1).  We  will  later  see  that  one  of  the  important  implementation 
challenges  is  the  task  of  efficiently  managing  the  structure  and  nonsingularity  of  i7^- 

The  full  factored  row  basis  is  then 


H 


n-(k+l) 


(J) 

En 

E\2 

F13 

}k 

(jj) 

F21 

F22 

F23 

}l 

{jjj) 

{  0 

0 

-I  ) 

}n- 

-(k  +  l) 

Introducing  the  notation 
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Au     =    E\  1  -  EnF22  F'2\ 
An    =    E\z  —  E\2F2~2  -^23 

where  A\\  is  the  Schur  complement,  or  Gauss  transform,  of  F22  in  A\\  (e.g.,  Golub  and  Van 
Loan  [1983]),  we  can  write  its  inverse  as 


B-'  = 


■Au  EuF22 


ATM 


11  ^13 


-F^F^ATl     (I  +  F^FvA^E^F^     Fnl{Fn-F2iAulAi3) 
0  0  -/ 


Grouping  the  coefficients  of  the  nonbasic  constraints  and  applying  the  same  column  ordering 
yields 


/  n-(k+l) 


D  = 


(i) 

-I 

0 

0 

}k 

(ii) 

0 

-/ 

0 

}l 

(Hi) 

F31 

F32 

£33 

}  m  —  k 

(iv) 

^F4] 

F42 

F43 

) 

}p-l 

The  principal  part  of  the  factored  tableau  is  DP,  where  P  =  —B    Ms  the  conjugate  row 
basis.  With  the  additional  notation 


131 


E 


31 


EwFnn    Fo 


-^33      =      £33  —  £32^22  Ex 


^32-r22  r21 
722ljP23 

F41     =    F41  -  F^F22  F21 
F43    =    F43  —  F42F2^  F23  , 


the  principal  part  of  the  factored  tableau  is 


DP  = 


(0 

(ii) 

(Hi) 
(iv) 


0)  Ui)  ijjj) 

i4n  -Ajj  Fi2F22  An  A13 

-F^Fajin3     (/ +  F221F21in1^i2)F221     F221(F23  -  F21iH1i13) 


[-1 


-A*\A[{         (Az\A[lEn  -  E32)F22l  A33  -  Asi^Mkj 


■Aiiifi1  (F4,iri1^i2-F42)F 


22 


F43-F4,A7,M 


41/ln  -4i3        / 


12 


Partitioning  w  =  (w),  w2,  io3,  u;4),     rT  =  (rj \r%)  and  bT  =  (bj,b%)  and  introducing  the 
notation 


W2 

= 

w2  -  vo\F22  F21 

W3 

= 

W3  -  W2F22F23 

h2 

= 

b2  —  F4XF22  b] 

f] 

= 

n-EuFn% 

h 

= 

r2  -  E2xF22]bx 

the  complete  factored  tableau  is 


1-1 


-  1 


■F22  F2XAU 
-isiifi1 

-F41AT1 

-w2Aux 


-An  E\2F22 

{I+F22'F2XA-xlEn)Fi2x 

{A^A^En-EZ2)F22x 

(F4XAU  E\2  -  F42)F22 

(w2AulEx2-wx)F22l 


AT?  Al3 


^22   (-^23 


-F2ii4nMi3) 

A3XA^AX3 


A33 

F43-FiXAnlAx3 


AXih 

F22\bx-F2XAuxfx) 
f2  -  A\\AX\f\ 

b2  -  F4XAjllr1 

h 


103  -  w2An  i4i3  wxF22bx  +  w2Aurx 


We  see  that  with  knowledge  of  the  current  factorization,  we  can  construct  the  entire  tableau 
from  F22  ,  A^  and  the  original  problem  data.  The  dimension  of  F22  is  equal  to  the  number  of 
F-type  constraints  that  are  currently  basic,  and  thus  can  be  at  most  p.  The  dimension  of  Axx 
is  equal  to  the  number  of  F-type  constraints  that  are  currently  basic,  and  thus  cannot  exceed 
m.  We  call  A^  the  explicit  transformation  kernel  and  F22  the  factored  transformation 
kernel. 

7     Factored  Column  and  Row  Generation 

Consider  generation  of  column  c  from  the  principal  part  of  the  tableau  DP.    Rewriting  in  a 
manner  that  highlights  our  intentions 


(7) 


Uj) 


(0 

DP  =  (i?) 
{Hi) 

(iv) 


[A] 


[-Au  EX2F22  } 


Uii) 

1 


[An1  An 


{-F22]F2X[  ]}        {F2-21-F221F21[  ]}     {F221F23-F221F21[  ]} 
-F3J  ]  -  F32{   }      — F3i[  ]  -  F32{   }      F33  -  F3][  ]  -  F32{    } 
V-F41[  ]-F42{   }       -F4][  ]-F42{    }       F43-F4i[  ]-F42{    }  ) 


where  Ax3  is  defined  as  before,  and  the  brackets  "[  ]"  and  "{  }"  contain  terms  common  to  but 
displayed  only  once  for  each  column. 

Assume  we  want  to  place  column  c  into  work  array  z.  We  partition  z  conformably  as 
zT  —  (zj^.z^z^zj),  refer  similarly  to  components  of  unit  vector  ec,  and  employ  a  (m  +  p)- 
vector  work  array  z.  The  notation  "<— "  denotes  simple  assignment,  "="  indicates  that  a  set  of 
factored  equations  must  be  solved,  and  "=*•"  explains  the  corresponding  result. 

If  column  c  is  in  (j), 
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Z\ 


K,'lc, 


and  then 


S<--F2i[zi] 

F22Z2  =  z 


Z2+-{-F221F21[z1}}, 


then 


and  finally 


if  c  is  in  (jj),  solve 


z3  < £31  [^i]  -  -^32(22}, 


z4*--Fu[z1]-F42{z2h 


F22Z2  =  ~(e2)c 

Z  <—  E\2Z2 


z  < (EpF22  )c 

z\  *-  [-Afj  EUF22  ]c, 


then 


z  <-  (e2)c  -  F2i[zi] 
F22z2  =  z 


^-{(F22l)c-F22lF21[zl}}, 


and  finally 


if  c  is  in  (jjj), 


Z3  < £31  [zi]  -  -£^32  {^2} 

Z4  < F41  [Z\]  -  -f42{-2}', 


*  «-  (^23)C 
F22Z2  =  Z 

z  <—  (Em)0  -  E12Z2 

Zi    4-    4_1- 


z2^F221(F23)c 
^-(il3)C 

zi  -  [irMiais 


then 


-  ♦—  (^23)°  -  F21Z1 
F22Z2  =  z 


^-{^'(^r-F^^.u-,]}, 


and  finally 
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~3  +—  (-^33)C  -  ^31  [*l]  -  E32{z2} 

Similarly,  row  r  can  be  generated.  The  principal  part  of  the  tableau  is  now  viewed  as 


DP  = 


0 

(«) 

(iii) 

(iv) 


/ 


0) 


Mtf] 


Ujj) 


\ 


[  ]J513  +  {    }F23 
h-P^Faiiin1]       {!#  -  [  \El3Fg)  [  ]£713  +  {   }F23 

[-iisiiiFl1]  {("^32  -[    l^)^}       £33  +  [    ]£l3  +  {    1^23 

V     [-AiiFi1]         {(-F42-[  JJSw)^1}     F43  +  [  ]£i3  +  {   }F2sJ 


where  ^31  and  F41  are  defined  as  before,  and  the  brackets  "[  ]"  and  "{  }"  contain  terms  common 
to  but  displayed  only  once  for  each  row. 

Assume  we  want  to  place  row  r  into  work  array  z.  We  partition  z  conformably  as  z  — 
(25,26,2:7),  refer  similarly  to  components  of  unit  vector  er,  and  employ  a  (m  +  p)-vector  work 
array  z . 


If  row  r  is  in  (i), 


l-^n  J') 


then 


z  < — hEn 
z%F22  =  z 


^^{-[An^rEuFv1}, 


and  finally 


27  <—  [^5]Fi3  +  {z6 }F2^; 


if  r  is  in  (ii),  solve 


•=6-^22  =  "Mr 
Z  <—  56-^21 

h  -  -zAul 


~{F22\ 
-(F22])rF2l   i 


then 


(e6)r  -  [h]En 


z&F22 


z6*-{(F22l)r-[zs)El2F221}, 


and  finally 
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if  r  is  in  (Hi), 


then 


and  finally 


if  r  is  in  (it>), 


then 


and  finally 


£7  «-  [h\E\z  +  {ie}^; 


2  «-  ~(^32)r 

^6-^22  =  2  =>       Z6  < (E32)rir22 

2  «-  (jRsi)r  +  *B*21        =>       2  <-  (Al)r 


5  < (i?32)r  -  [25] -^12 

56F22  =  z  =»     f6  -  {(-(£32)r  -  N^)^1} 


*r  <-  (^33 )r  +  [5s]Ei3  +  {i6}^23; 


S<--(*42)r 

^6-^22  =  -  =*•        h  < -(-f42)r^22' 

Z  *-  (F4l_)r  -  i6*21        =>        Z  *-  (Al)r 


-  < (-^42)r  -  [-5] -El 2 

Z6F22  =  Z  =*        =6  -  {(-(E42)r  -  [^fi^)^'}' 


57^(i?43)r  +  [i5]Ei3+{£6}E: 
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8      The  Complete  Algorithm 

The  complete  algorithm  is  described  in  terms  of  abstract  functions  which  operate  on  funda- 
mental data  structures. 

Tableau  management  requires  two  index  maps:  one  yielding  the  intrinsic  coordinate  in  the 
principal  part  of  the  tableau  for  each  original,  extrinsic  problem  row  or  column,  and  the  other 
its  inverse  map.  Intrinsic  arguments  are  shown  in  lower  case,  and  extrinsic  in  upper  case. 
Index _Exchange(indexl,index2)  updates  these  maps  for  the  exchange  of  a  pair  of  tableau 
coordinates  indexl  and  index2. 

The  tableau  regions  are  successive  partitions  of  indices 
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TABLEAU 

ENDING 

REGION 

INDEX 

(») 

MEC 

(«) 

MFC 

(m) 

MER 

(w) 

r// 

0) 

NER 

0j) 

NFR 

(jjj) 

m  +  n 

CONTENTS  OF  REGION 

basic  Columns  solving  Explicit  rows 

basic  Columns  solving  Factored  rows 

nonbasic  Explicit  Rows 

nonbasic  Factored  Rows 

basic  Explicit  Rows 

basic  Factored  Rows 

nonbasic  Columns. 

Increment (endingJndex)  and  Decrement(endingJndex)  are  functions  to  modify  these 
ending  indices. 

Generate_Row(row)  and  Generate_Column(column)  place  numeric  values  of  a  tableau 
row  or  column  in  ROWCOL0,  which  is  commensurate  with  the  tableau  dimensions. 

Using  ROWCOL(),  Update_Rim(row,col)  maintains  current  numeric  values  of  the  right- 
hand-side  and  bottom  row  of  the  tableau  in  RIM(). 

The  explicit  transformation  kernel,  A^i  ,  is  operated  on  by  functions  using  ROWCOLQ: 

Add_E-Kernel_Row(ROW), 
Add_E-Kernel_Column(COLUMN) , 
Delete_E-Kernel_Row(ROW), 
Delete_E-Kernel_Column(COLUMN), 

Replace_E-Kernel_Row(REPLACED_ROW, REPLACING _ROW),  and 
Update_Explicit _Transforrnation_Kernel,  the  pivotal  update. 

A^  can  be  represented  in  any  way  that  suits  the  implementer  and  efficiently  supports  these 
functions.  Generally,  we  find  .AT,1  to  be  relatively  sparse,  and  report  here  results  using  an  array 
with  an  entry-point  for  each  inverse  row  and  column  which  accesses  a  stack  of  orthogonally- 
linked  nonzero  inverse  elements.  We  have  also  implemented  a  dense  vector-processor  version, 
and  the  design  issues  for  LU-based  schemes  are  given  by  Olson  [1989]. 

The  factored  kernel,  F22,  is  manipulated  with: 

Factored_Kernel_Singular(ROW, COLUMN),  a  Boolean  function, 
Find_E-Kernel_Column_for_Key(ROW,COLUMN),  a  column  from  region  (i), 
Find_F-Kernel_Column_to_Remove(ROW, COLUMN),  a  column  from  region  (ii),  and 
Update_Factored_Kernel(ROW, COLUMN),  the  pivotal  update. 
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The  complete  abstract  algorithm  is: 


0.  Initialize; 

1.  select  primal  or  dual  algorithm; 

2.  Select  a  primal  (col)  or  dual  (row)  violation; 
STOP  if  current  solution  is  terminal,  or 
Generate_Column(col)  or  Generate_Row(row); 

3.  By  ratio  test  with  RIM()  and  ROWCOL(), 
select  a  (row)  or  (col)  pivot  coordinate; 
STOP  if  current  solution  is  terminal,  or 
Generate_Row(row)  or  Generate_Column(col); 

4.  Update_Rim(row,col), 

primary  Index_Exchange(row,col) , 

perform  secondary  and  tertiary  exchanges  (Table  8-1), 

Update_Factored_Kernel(ROW,COL) 

Update_Explicit_Transformation_Kernel, 

Go  to  step  1. 


18 


;=;< 


s- 
o 
•at 

-i 
O 
O 


•       3 

a    o 
a    * 

I    £      ~ 

4    o    „« 

i^n 2  2s 

J  "5  >--^-'<=-f--- 

s-'  q  9  9  «  y 

,    E    M   M    k± 

-  a  - 


L?H 


i    3 


£■ 

o 
£? 

4* 

Ml 

isgl 

<i  84 

LI  111 


2? 


■g 


o  a 

o  - 


il 

»  a 

I  3 

1  o 


ST 
0 

z 

Ti   2 


1       £z-2 

!  s'z  * 

1  ^     •£.&  o 

!  a:  a?  •  •  "J 

!  CiJ  u     !f  !f  — 
1  Z  2    S    9     * 

1  C-O  !  _•   a 

!    «•    «•  *g    t»     « 

'  a  a  x  m  ^ 

1     •     9  fj  ul 


u  u  •  « 

S5.  S  3 


*   «   • 
C  -o  -a 


5 
o 

o 
u 


—  !    0    O    Q  ■■  ■» 


«    S'   Sa.«   •   9    S   G  * 


I  °4 

a  „  -o 
2,2  5 


•  4  (i  u 

id    •  4    * 

5  •  •  5  S 

£os:.2g 

u 


ssss 

J  -c  _  _ 
«  .  fl  ° 

44  s  1 


o 


4  _ 

5  So 
1  oSJ 

3  *■* 


.2       3  ■2'?  * 


4 

■s-s  g 


Oa? 

HI 

^  e  x 

4§3 

<.35 


—     3      *     S  « 

s  fl-.  ■■  -fl     -a  — 

^  »  v  «   a   k   : 

4i^n 

'  a  -a  S  S  •  ■£ 


2 


»1 


Hi 


^■3  - 
i  n  §ii  84 


■*  a: 

IS 


1  '*"  s 
1  wi*  a:  -  2  5 


!^ 


2f 

3J  O 

oB 

J* 

-J   0 

o  a 

O    v 

X«l 

la 

H  a 

■s  i- 


4s 

12 


£s 

q   o  —  O 

■iri    • 

«  •  a 
-1  -  .  *  _ 


a   a 


*  1 


A  =  7  5 
*i!ti  u 


»  as  •  ss   (  « 
o  "O  a  "9  a.  •  fe 


u. 


Id 


a 

03 

u 
X 

U 

3 

03 

6 


1- 

-a 
c 

03 

03 

-a 
c 
o 
u 
1) 


3< 


iMl 


■o-o  i 


a  *-f  f 

*  S  °  3  a  s  s 

—  a   i,   «   «  —  — 

-  :  v  -o  j  -  - 

llS^44li 

*      4'       1     —     *     X     *     * 

^  St  •  •  •  t  5 

*  3  9  "2  "2  5  * 

2u3i:oa 


0 

Bos 


31 

Mi 


g  a  z 

air 

^4  i 

-  *  £ 
•  1  8 


O  a?" 


a  -1   "=  1 

i s ISg s 

A  x   x    J   a   k 

-  «*  *3  i  1^ 

J  j  k  t  ;  « 

*   *C     ■     g     Z     9 

»  J"s  *  a "g 


00 

W 
-J 
ffl 


u    a> 


.a  s 


19 


Functions  for  maintaining  the  factored  kernel  vary  with  the  factorization,  and  a  good  im- 
plementation will  exploit  these  differences.  However,  a  general  specification  will  suffice  for  all 
factorizations  here. 

We  are  exclusively  interested  in  performing  two  fundamental  operations: 

1.  Solving  factored  equations  of  the  form: 

F22Z2  —  t>2 

and 

^2  F22  =  b2 
where  Z2  and  Z2  are  unknown  and  62  and  62  are  rational  (not  necessarily  integer),  and 

2.  Restoring  F22  to  the  desired  form  of  F22  which  makes  the  factored  equations  above  easy 
to  solve,  where  F22  results  from  inflicting  F22  with  a  column  exchange,  a  column  and  row 
deletion,  a  column  and  row  addition,  or  a  row  exchange. 

Factored_Kernel_Singular(LABELl,LABEL2)  predicts  whether  F22  will  be  singular  if: 

1.  LABEL1  gives  a  basic  column  in  a  basic  factored  ROW  for  which  COL=LABEL2  would 
be  exchanged; 

2.  LABEL1  is  a  basic  factored  ROW  solved  by  basic  column  COL=LABEL2,  both  of  which 
would  be  removed  from  the  factored  kernel; 

3.  LABEL1  is  a  nonbasic  factored  ROW,  which  would  be  added  to  the  factored  kernel  with 
column  COL=LABEL2,  or 

4.  LABEL1  is  a  basic  factored  ROW  which  would  be  replaced  by  some  other  row  LABEL2. 

The  first  case,  a  column  exchange,  is  equivalent  to  asking  whether  solving  F22Z2  —  ^2  with  62 
equal  to  region  (ii)  of  the  proposed  entering  column  COL,  yields  a  nonzero  term  associated  with 
the  basic  factored  row  in  Z2(row).  This  is  easy  to  answer  for  the  factorizations  we  discuss.  For 
instance,  Brown  and  McBride  [1984]  show  that  back-solving  a  "nearly-triangulated  generalized 
network"  basis  F22  only  as  far  as  ROW  will  suffice,  and  that  this  can  be  implemented  as  a 
traversal  of  the  "backpaths"  of  the  (zero,  one,  or  two)  coefficients  in  62  for  the  column  now 
basic  in  ROW.  If  numerical  precision  is  an  issue,  Z2(row)  can  be  computed  during  this  search 
and  tested  for  significance. 

The  second  case,  a  row  and  column  deletion,  is  trivial  because  the  resulting  F22  will  always 
be  nonsingular. 

The  third  case,  a  row  and  column  addition,  can  be  answered  by  adding  ROW  and  COL  to 
F22  without  disturbing  its  desirable  special  structure,  thus  creating  an  instance  of  case  one.  For 
instance,  placing  ROW  first,  and  COL  last  in  F22  suffices  for  all  factorizations  discussed  here. 

The  fourth  case,  a  row  exchange,  can  be  answered  by  applying  the  second  case,  deleting 
ROW  and  its  basic  column,  and  then  (perhaps  repeatedly)  applying  the  third  case,  adding  the 
row  with  index  COL  and  any  column  which  would  yield  a  nonsingular  F22- 

Find_F-Kernel_Column_to_Remove(ROW,COL)  given  basic  factored  row  ROW,  re- 
turns its  basic,  or  "key"  column  COL. 

Find_E-Kernel_Column_for_Key(L ABEL, COL),  given  either  a  nonbasic  factored  ROW= 
LABEL,  or  a  basic  column  LABEL  solving  a  basic  factored  ROW,  searches  for  an  acceptable  ba- 
sic column  COL  in  Explicit  Rows  (region  (i))  using  Factored_Kernel_Singular(ROW,COL). 

Update_Factored_Kernel(LABELl,LABEL2)  restores  F22  to  the  desired  form  of  F22, 
with  a  possible  increase  or  decrease  in  dimension.  Following  the  case- by-case  scheme  of  Fac- 
tored_Kernel_Singular,  we  pre-  and  post-process  the  factored  basis  representation  to  permit 
use  of  a  single,  static  factorization  update  of  conventional  design.  Olson  [1989]  pursues  this  in 
considerable  detail. 

Table  8-2  displays  supporting  data  structures  for  these  functions. 
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TABLE  8-2  Factorization  Algorithm  Data  Structures 

DATA 
FACTORIZATIONS     STRUCTURE      SIZE      USE 


GUB,PN,GN 

RIMQ 

m  +  n 

ROWCOL() 
MSKRC() 

m  +  n 
m  +  n 

LQRC() 

m  +  n 

KEY() 

m  +  n 

PN,GN 

WORK() 

P 

MSKWK() 

V 

LQWK() 

V 

PO() 

P() 

V 
V 

D() 


current  tableau  right-hand  side, 
bottom  row 

current  tableau  row  and  column 
logical  mask  true  for  corresponding 
nonzero  in  ROWCOL(),  false  otherwise 
LIFO  queues  of  nonzero  row  and 
column  coordinates  in  ROWCOLQ 
basic  column  in  basic  factored  row, 
and  vice  versa 

values  for  62  or  62 

in  basic  factored  equations 

logical  mask  for  nonzero  row  and 

column  coordinates  in  WORK() 

LIFO  queue  of  nonzero  coordinates 

in  WORK() 

next  basic  factored  row  in  pre-order 

off-diagonal  row  with 

nonzero  factored  coefficient 

depth,  remaining  back-substitution 

path  length  in  factored  component 


GN 


VGN() 
JMULQ 


p         generalized  network  cycle  factors 
n  ratio  of  generalized  network 

coefficients  in  each  column 


Generalized  upper  bound  (GUB),  pure  network  (PN),  and  generalized  network 
(GN)  factorizations  respectively  require  more  data  structures  to  support  kernel 
factorization.  For  each  factorization,  F22  is  maintained  in  some  partial  ordering 
of  rows,  and  of  columns  —  a  signed  identity  for  GUB,  upper-triangular  for  PN 
(e.g.,  Bradley,Brown,  and  Graves  [1977]),  and  nearly-upper-triangular  for  GN  (e.g., 
Brown  and  McBride  [1984]).  Direct  solution  of  factored  kernel  equations  F22Z2  =  ^2 
and  z%F22  —  by  ls  performed  alternately  using  the  data  structures  shown. 
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9      Computational  Experience 

The  factorization  methods  introduced  here  have  been  implemented  and  used  to  solve  a  variety 
of  existing  models  provided  by  our  colleagues.  With  their  help,  we  have  extracted  suggested 
problem  instances  from  a  diversity  of  decision  support  systems  on  host  computers  ranging 
from  mainframes  to  micro-computers.  Because  our  goal  is  to  test  factorization  technology  in 
isolation,  the  results  reported  here  are  achieved  without  benefit  of  any  model-specific  knowledge 
or  tuning. 

However,  we  also  seek  to  develop  effective  modeling  tools  for  customized  use  in  developing 
and  refining  new  models.  To  this  end,  we  have  greatly  benefited  from  the  experience  and  advice 
of  our  colleagues,  and  we  include  some  discussion  of  modeler  guidance  and  insight  along  with 
the  numerical  results. 

Each  model  is  introduced  below  by  a  short  synopsis.  Multiple  instances  of  some  models 
are  reported  where  diversity  of  size,  structure,  and  taxonomy  have  proven  interesting  to  the 
modelers.  For  those  models  that  employ  nonlinear,  mixed  integer,  or  decomposition  features, 
we  report  solution  statistics  for  the  initial  linear  program. 

•  GTE  The  seven  Telephone  Operating  Companies  within  GTE  have  adopted  an  integrated  business  system 
called  Capital  Program  Management  System  (CPMS)  to  guide  their  3  billion  dollar  per  year  capital  planning. 
The  system  includes  a  large  scale  mixed  integer  programming  optimization  system  that  optimizes  the  critical 
economic  tradeoffs  between  maximizing  the  long-term  budget  value  of  the  firm's  equity  and  satisfying  shorter- 
term  financial  constraints,  resource  limitations  and  service  objectives.  Investment  opportunities  for  the  next 
5  years  are  modeled  as  0-1  variables  with  alternative  implementations  for  each.  The  objective  is  to  maximize 
the  net  present  value  of  the  capital  portfolio.  There  are  financial  constraints  on  capital,  internally  generated 
funds,  net  income  to  common,  and  limits  on  resources  such  as  labor  hours,  lines  installed,  etc.  There  are  also 
constraints  that  enforce  logical  relationships  among  opportunities  (such  as,  if  choose  A  then  must  choose  B). 
See  Bradley  [1986]. 

•  INVEST  Capital  allocation  and  project  selection  for  Mobil  Oil  Corporation  are  modeled  as  a  two-stage 
multi-year  nonlinear  capital  budgeting  problem  with  over  40,000  integer  variables.  A  master  problem  allo- 
cates capital  among  markets  over  a  multi-year  horizon  considering  the  estimated  nonlinear  effects  on  sales 
of  concentrated  marketing  investments.  The  instance  reported  here  is  a  mixed-integer  linear  program  sub- 
problem  of  the  two-stage  model  which,  given  these  annual  capital  expenditure  limits  for  a  market,  selects 
particular  alternate  investments.  Such  subproblems  are  easy  to  solve,  and  optimality  is  achieved  with  a  single 
iteration  of  the  nonlinear  master  problem.  See  Harrison,  Bradley  and  Brown  [1989]. 

•  TANKER  A  crude  oil  tanker  scheduling  problem  faced  by  a  major  oil  company  is  solved  using  an  elastic 
set  partitioning  model.  The  model  takes  into  account  all  fleet  cost  components,  including  the  opportunity 
cost  of  ship  time,  port  and  canal  charges,  demurrage  and  bunker  fuel.  The  model  determines  optimal  speeds 
for  the  ships  and  the  best  routing  of  ballast  (empty)  legs,  as  well  as  which  cargoes  to  load  on  controlled  ships 
and  which  to  spot  charter.  All  feasible  schedules  are  generated,  the  cost  of  each  is  accurately  determined 
and  the  best  set  of  schedules  is  selected.  See  Brown,  Graves  and  Ronen  [1987). 

•  HFDF  A  large-scale  elastic  set  partitioning  model  used  to  assign  frequencies  for  a  network  of  high  frequency 
direction  finding  receivers.  See  Brown,  Drake,  Marsh,  and  Washburn  [1990]. 

•  GAS  A  multi-time  period  strategic  model  for  use  by  natural  gas  utilities  for  determining  optimal  contract 
levels  for  gas  purchase,  storage  and  transmission.  An  underlying  generalized  network  flow  model  represents 
gas  being  bought,  stored,  shipped  and  consumed  over  a  multi-year  time  horizon,  typically  at  a  monthly  level 
of  detail.  Constraints  and  variables  are  added  to  handle  variable  maximum  and  minimum  purchase  levels, 
variable  leased  or  constructed  storage  and  variable  transmission  capacities.  An  integrated  parallel  model 
incorporates  the  peak  requirements  necessary  on  some  days  during  cold  winter  months.  This  model  has  been 
used  by  a  number  of  utilities  including  Southwest  Gas  Corporation  and  Questar  Pipeline  Corporation  to  plan 
operations  and  to  justify  such  plans  to  regulatory  agencies.  See  Avery,  Brown,  Rosenkranz  and  Wood  [1992]. 

•  KELLOGG  A  multi-time  period,  multi-plant  production/inventory/transshipment  linear  program  for  Kellogg 
cereals.  The  model  guides  weekly  processing,  packaging  and  shipping  decisions.  Production  consists  of  two 
stages:  processing  lines  produce  basic  products  which  are  then  packaged  on  packaging  lines  into  different- 
sized  containers  to  yield  finished  products.  Processing  lines  produce  a  subset  of  the  basic  products  and  have 
limited  capacity  with  overtime  charges  for  weekend  shifts.  Packaging  lines  are  analogous.  In-house  inventory 
capacity  is  limited  although  outside  storage  is  available  at  additional  cost.  Inter-plant  shipments  of  finished 
products  are  made  by  rail  or  truck.  See  Wood  [1989]. 
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•  ODS  A  commonly  occurring  problem  in  distribution  system  design  is  the  optimal  location  of  intermediate 
distribution  facilities  between  plants  and  customers.  A  multicommodity  capacitated  single-period  version 
of  this  problem  is  formulated  as  a  mixed  integer  linear  program.  A  solution  technique  based  on  Benders 
Decomposition  is  developed.  ...  An  essentially  optimal  solution  is  found  and  proven  with  a  surprisingly  small 
number  of  Benders  cuts.  See  Geoffrion  and  Graves  [1974].  The  instances  reported  here  are  decomposition 
master  problems. 

•  TAM  The  annual  decision  on  how  much  the  Air  Force  should  spend  on  aircraft  and  on  munitions  is  of  great 
interest  to  many  people.  How  the  Air  Force  staff  develops  information  to  support  the  decision  has  changed 
over  the  years.  Currently,  a  linear  program  is  being  used  by  the  Air  Force  Center  for  Studies  and  Analysis 
and  is  being  tested  by  the  Munitions  Division  of  the  Plans  and  Operations  Directorate  (AF/XOXFM)  for 
municions  tradeoff  analysis.  The  LP  uses  existing  data  and  estimates  of  (1)  aircraft  and  munition  effectiveness, 
(2)  target  value,  (3)  attrition,  (4)  aircraft  and  munition  costs,  and  (5)  existing  inventories  of  aircraft  and 
munitions.  Other  factors  considered  are  weather  and  length  of  the  conflict.  See  Might  [1987]  and  Jackson 
[1989]. 

•  PHOENIX  A  planning  model  for  the  multi-year,  multi-billion  dollar  modernization  of  the  U.  S.  Army's  aging 
helicopter  fleet.  The  mixed  integer  linear  program  employs  a  multi-product  production/inventory  formulation 
with  aged  inventory.  Goal  constraints  attempt  to  enforce  fleet  size,  maximum  age,  and  technology  goals  for 
each  year  and  each  of  four  aircraft  missions,  while  also  keeping  expenditures  within  upper  and  lower  limits. 
Additionally,  combinatorial  constraints  and  variables  handle  production  line  startup  and  shutdown  costs, 
minimum  and  maximum  production  levels  and  requirements  linking  certain  production  lines.  See  Clemence, 
Teufert,  Brown  and  Wood  [1988]  and  Brown,  Clemence,  Teufert,  and  Wood  [1990]. 

•  EA6B  Configures  jammers  of  hostile  radar  on  an  EA-6B  "Prowler"  Naval  electronic  warfare  aircraft.  See 
Sterling  [1990]. 

•  DEC  Digital  Equipment  Corporation  uses  this  model  to  determine  worldwide  manufacturing  and  distribution 
strategy  for  new  products.  This  mixed  integer,  linear  program  suggests  a  production,  distribution,  and  vendor 
network  which  minimizes  cost  and/or  cumulative  cycle  times  subject  to  constraints  on  estimated  demand, 
local  content,  and  joint  capacity,  over  multiple  products,  echelons,  and  time  periods.  Cost  factors  include 
fixed  and  variable  production  charges,  distribution  via  multiple  modes,  taxes,  duties  and  duty  drawback,  and 
inventory  charges.  See  Harrison,  Arntzen,  and  Brown  [1992]. 

•  AMMO  4H  A  four-commodity  transshipment  model  for  delivery  over  time  of  military  products  from  pro- 
duction and  storage  locations  to  overseas  locations  to  support  theater  operations  is  developed.  The  model 
covers  five  physical  echelons,  including  production  plants,  storage  depots,  ports  of  embarkation,  ports  of 
debarkation  and  geographic  field  locations.  Road,  rail,  sea  and  air  transportation  are  modeled,  and  product 
demands  are  time-phased.  Capacitation  occurs  primarily  on  sea  and  air  links,  and  as  throughput  capacities 
on  transfer  points,  requiring  replication  of  some  echelons.  The  objective  of  the  model  is  to  minimize  deviation 
from  on-time  deliveries.  See  Staniec  [1984]. 

•  BUSCH  A  model  of  brewery-to-wholesaler  movements  of  beer  for  Anheuser-Busch.  The  model  also  includes 
some  packaging  decisions  and  is  essentially  a  multicommodity  flow  model  with  joint  capacity  constraints 
arising  from  loading  dock  and  inventory  capacities  as  well  as  some  managerial  requirements.  See  Brown, 
Mamer,  McBride  and  Wood  [1992].  The  instance  reported  here  is  a  small  pilot  model  for  the  full-scale 
system  with  millions  of  variables  which  is  solved  directly,  or  by  decomposition. 

•  BAR  A  linear,  mixed-integer  multi-period  production-inventory  master  planning  model.  See  Harrison  [1992]. 

Four  implementations  are  compared:  "XS"  is  an  unadorned  version  of  the  X-System,  an 
implementation  of  the  Graves  mutual  primal-dual  method  with  its  GUB  factorization  disabled, 
while  "XS(GUB)",  "XS(PN)",  and  "XS(GN)"  each  employ  the  respective  factorizations  dis- 
cussed here.  To  estlablish  a  frame  of  reference,  performance  of  these  implementations  is  com- 
pared with  two  well-known  commercial  solvers:  IBM's  Optimization  Subroutine  Library  "OSL" 
(Release  2  [1991]),  and  two  versions  of  "CPLEX"  (Version  1.2  [1990]  and  Version  2.0  [1992]). 

Ideally,  one  would  develop  four  equivalent  formulations  of  each  model,  each  customized  for 
its  particular  solver  with  the  goal  of  inducing  a  large  factored  row  set  of  the  appropriate  type. 
This  approach  is  a  consistent  theme  in  the  literature  dealing  with  specialized  algorithms  and 
one  that  we  strongly  endorse.  Alternate  formulations  of  a  model  are  often  available,  and  it 
seems  sensible  to  choose  one  that  exploits  as  much  as  possible  the  strengths  of  the  solver. 
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However,  all  of  the  models  used  here  are  "off-the-shelf"  in  the  sense  that  they  were  developed 
at  various  times  by  various  modelers,  and  alternate  formulations  are  impracticable.  Thus,  the 
approach  is  to  preserve  a  single,  unfactored  representation  of  each  model,  and  attempt  to 
identify  favorable  row  structures  through  the  use  of  heuristics.  The  procedure  is  based  on 
the  work  of  Brown  and  Thomen  [1980],  Brown  and  Wright  [1983]  and  Brown,  McBride  and 
Wood  [1985].  The  heuristics  are  greedy  and  myopic  in  the  sense  that  they  initially  consider  the 
entire  row  set  of  the  problem,  and  discard  one  row  at  a  time  without  backtracking  until  the 
remaining  set  satisfies  the  desired  row  factorization.  This  can  be  expected  to  confound  or  destroy 
structure  introduced  by  the  modeler.  Although  the  automatic  factorization  implementation  has 
options  to  accept  modeler  guidance,  the  methods  are  compared  here  without  this  subjective 
complication.  While  these  model-naive  experiments  yield  interesting  and  useful  observations 
about  the  implementations,  they  suffer  for  lack  of  guidance  by  a  skilled  modeler. 
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Table  9-1  shows  the  important  structural  information  concerning  the  model  instances  to  be 
solved. 

TABLE  9-1  Problem  Dimensions 


n 

m 

Pgub/% 

Ppn/% 

Pgn/% 

NZEL 

GTE 

6,624 

960 

909/95 

909/95 

922/96 

58 

INVEST 

11,989 

1,338 

941/70 

1,101/82 

1,168/87 

33 

TANKER 

7,598 

83 

32/39 

32/39 

66/80 

31 

HFDF 

10,548 

61 

31/51 

31/51 

32/52 

189 

GAS  PN  A 

27,884 

6,848 

4,345/63 

5,934/87 

5,976/87 

37 

GAS  PN  C 

15,362 

3,794 

2,658/70 

3,418/90 

3,420/90 

20 

GAS  PN  E 

5,102 

1,184 

434/37 

877/74 

883/75 

7 

GAS  GN  A 

27,884 

6,848 

4,484/65 

5,142/75 

5,976/87 

37 

GAS  GN  C 

15,362 

3,794 

2,664/70 

3,084/81 

3,420/90 

20 

KELLOGG  2 

17,841 

3,818 

1,265/33 

2,578/68 

2,596/68 

35 

KELLOGG  3 

27,490 

5,727 

2,295/40 

3,867/68 

3,892/68 

54 

KELLOGG  4 

37,139 

7,636 

2,428/32 

5,156/68 

5,188/68 

71 

KELLOGG  5 

46,788 

9,545 

3,388/36 

6,445/68 

6,484/68 

93 

ODS  1 

11,568 

3,023 

528/17 

540/18 

558/18 

21 

ODS  3 

23,993 

594 

490/82 

490/82 

490/82 

68 

TAM  5 

10,531 

438 

102/23 

132/30 

162/37 

94 

TAM  8 

6,104 

420 

118/28 

154/37 

196/47 

49 

TAM  12 

17,793 

629 

177/28 

231/37 

294/47 

165 

PHOENIX  10 

6,884 

1,618 

206/13 

220/14 

1,153/71 

14 

PHOENIX  30 

17,212 

4,305 

293/07 

303/07 

3,604/84 

48 

EA6B 

12,247 

2,978 

1,610/54 

2,921/98 

2,921/98 

16 

DEC 

14,518 

2,171 

677/31 

677/31 

1,088/50 

24 

AMMO  4H 

83,497 

13,963 

6,874/49 

12,892/92 

12,892/92 

129 

BUSCH  4 

7,997 

1,248 

649/52 

1,140/91 

1,148/92 

15 

BAR 

49,032 

7,446 

2,712/36 

4,575/61 

5,134/69 

102 

For  each  problem,  the  total  number  of  structural  variables  is  n,  structural  con- 
straints m,  GUB  rows  found  by  the  identification  heuristic  pcan,  pure  network 
rows  ppn,  generalized  network  rows  Pcn,  aQd  thousands  of  nonzero  technological 
coefficients  NZEL.  For  example,  GTE  has  58  thousand  nonzero  technological  co- 
efficients for  6, 624  variables;  GTE  can  be  viewed  as  having  960  explicit  and  no 
factored  rows,  or  as  909/95%  GUB-factored  and  960  -  909  =  51  explicit  rows,  as 
909  PN-factored  rows,  or  as  922  GN-factored  and  38  explicit  rows. 
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Table  9-2  displays  solution  times  for  the  sample  problems.  These  CPU  times  exclude  initial 
problem  input,  factored  row  identification,  and  final  output  —  on  average  about  0.2  second  per 
problem. 

TABLE  9-2  Solution  Seconds 


GTE 
INVEST 
TANKER 
HFDF 
GAS  PN  A 
GAS  PN  C 
GAS  PN  E 
GAS  GN  A 
GAS  GN  C 
KELLOGG  2 
KELLOGG  3 
KELLOGG  4 
KELLOGG  5 
ODS  1 
ODS  3 
TAM  5 
TAM  8 
TAM  12 
PHOENIX  10 
PHOENIX  30 
EA6B 
DEC 

AMMO  4H 
BUSCH  4 
BAR 


AMDAHL  5995-700 

486/33MH2 

:  PC 

X-Syst 
GUB 

em 
PN 

GN 

IBM 

OSL 

XS 
GN 

CPLEX 

none 

1.2 

2.0 

9 

8 

8 

8 

43 

41 

65 

58 

4 

5 

4 

4 

23 

24 

52 

49 

8 

10 

10 

8 

6 

48 

8 

9 

100 

101 

101 

111 

40 

462 

581 

542 

2,376 

778 

50 

57 

291 

321 

1,714 

1,221 

4 

4 

3 

3 

79 

26 

185 

245 

72 

33 

4 

6 

13 

33 

165 

50 

1,115 

542 

294 

72 

312 

385 

3,532 

2,262 

4 

4 

3 

3 

71 

26 

186 

236 

3 

2 

2 

2 

69 

5 

219 

194 

5 

5 

5 

5 

144 

9 

513 

440 

48 

36 

26 

24 

500 

115 

1,274 

1,111 

1,122 

1,459 

320 

248 

1,210 

926 

2,615 

2,243 

6 

3 

2 

5 

125 

22 

1,042 

414 

6 

6 

6 

6 

23 

38 

58 

56 

55 

60 

45 

40 

44 

231 

151 

86 

24 

14 

14 

17 

21 

69 

77 

71 

108 

112 

88 

106 

101 

312 

416 

378 

4 

2 

2 

2 

20 

10 

84 

73 

41 

25 

26 

9 

335 

69 

1,171 

1,239 

75 

33 

9 

9 

45 

44 

177 

244 

189 

32 

42 

80 

71 

93 

317 

293 

42 

41 

35 

46 

1,000 

277 

1,762 

1,597 

5 

5 

4 

4 

26 

34 

55 

46 

290 

212 

106 

114 

382 

480 

1,366 

1,222 

An  AMDAHL  5995-700  running  under  IBM  VM/CMS/XA  with  IBM  VS  FOR- 
TRAN 2.3.0  is  used  to  render  performance  in  CPU-seconds  accurate  to  the  precision 
shown  for  the  basic  "XS"-system  using  no  factorization,  compared  with  dynamic 
factorizations  of  "XS(GUB)",  pure  network  "XS(PN)",  and  generalized  network 
"XS(GN)"  rows.  "IBM-OSL"  shows  primal  simplex  performance  on  the  same  com- 
puter of  the  IBM  Optimization  Subroutine  Library,  Release  2  [1991].  "486/33" 
shows  the  clock-time  performance  of  "XS(GN)"  on  a  microcomputer  (33  MHz  Intel 
486  with  32  MB  RAM)  followed  by  that  of  "CPLEX"  (Version  1.2)  [1990]  and  of 
"CPLEX"  (Version  2.0)  [1992]  on  the  same  microcomputer. 


The  original  formulations  of  most  of  the  test  problems  are  strongly  influenced  by  the  mod- 
eler's solution  strategy.  For  instance,  TANKER  is  endowed  with  a  GUB  structure  which  places 
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every  binary  variable  in  an  associated  set  from  which  only  one  member  can  be  chosen;  this 
maximal  GUB  set  is  also  sought  for  its  tendency  to  yield  nearly-integer  solutions  to  the  linear 
program.  AMMO  4H  is  a  multicommodity  capacitated  transshipment  problem  and  thus  is  best 
suited  to  a  pure  network  factorization;  it  was  originally  solved  by  dual  decomposition  rendering 
pure  network  sub-problems.  PHOENIX  is  a  multicommodity  equipment  replacement  model 
closely  following  the  generalized  network  factorization  paradigm. 

The  row  structures  in  Table  9-2  have  been  found  without  modeler  help  by  our  automatic 
identification  heuristic,  yet  they  corroborate  the  modelers'  intentions.  TANKER  reveals  the 
same  pure  network  rows  as  the  GUB  rows,  or  a  generalized  network  constructed  by  identifying 
one  additional  row  to  be  paired  with  each  GUB  row.  PHOENIX  exhibits  a  dominant  structure 
that  is  clearly  a  generalized  network. 

One  would  anticipate  the  factorization  exploiting  the  dominant  row  structure  to  win  com- 
putation tests.  This  is  wrong  more  often  than  right.  Table  9-2  shows  that  the  more  general 
factorization  dominates  the  less  general,  with  few  exceptions:  GTE's  relatively  large  GUB  set, 
and  the  large  pure  networks  in  GAS  PN  A,  GAS  PN  E,  and  AMMO  4H  seem  to  satisfy  our 
prior  bias  toward  model-dictated  factorizations. 

Table  9-2  also  suggests  that  our  myopic  use  of  heuristics  to  automatically  identify  factored 
structure  has  its  pitfalls.  In  a  number  of  problem  instances,  we  identify  significantly  larger  fac- 
tored sets  with  the  more  general  factorizations,  yet  we  enjoy  little  improvement  in  computation 
times  (e.g.,  INVEST,  ODS  3,  TAM  8).  This  suggests  that  the  "quality"  of  a  row  factorization 
is  not  completely  specified  by  its  size. 

We  have  pursued  this  notion  by  inviting  some  of  the  modelers  to  guide  our  identification 
heuristic  to  precisely  the  row  sets  they  intended.  Some  of  the  results  have  been  striking:  Wood 
[1989]  reports  significant  improvements  for  problems  in  the  GAS  system. 

A  few  of  our  corresponding  modelers  have  had  the  opportunity  to  build  models  from  scratch 
with  a  particular  factorization  in  mind.  This  admits  model  coercion  and  a  wide  range  of  well- 
known  reformulation  methods  which  we  think  can  materially  change  both  the  size  and  quality 
of  the  result.  Their  early  reports  show  promise.  Among  the  models  discussed  here,  the  GAS 
and  KELLOGG  systems  have  been  subsequently  reformulated  to  enhance  generalized  networks, 
GTE  has  been  re-engineered  to  further  accentuate  its  dominant  GUB  set,  and  TANKER-like 
and  many  ODS  models  have  been  moved  to  a  micro-computer;  all  these  models  are  now  larger, 
but  much  easier  to  solve. 

It  is  surprising  and  encouraging  that  the  transition  to  more  general  factorizations  seldom 
degrades  performance  much,  even  when  few  additional  factored  rows  are  won  by  the  increased 
generality.  This  contradicts  popular  folklore  that  the  more  general  factorizations  demand  sub- 
stantial, if  not  overwhelming  increases  in  the  resulting  sizes  of  the  factored  structures.  In 
fact,  computational  testing  reported  by  others  has  usually  been  limited  to  models  in  which 
the  number  of  explicit  rows  is  in  the  range  of  one  to  twenty  (e.g.,  Chen  and  Saigal  [1977], 
Glover,  Karney,  Klingman  and  Russell  [1978],  Glover  and  Klingman  [1981]).  Our  results  are 
all  the  more  remarkable  given  the  lack  of  guidance  from  the  modeler  for  the  "intended"  row 
factorization. 
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Table  9-3  shows  the  maximum  size  of  the  An   in  terms  of  its  nonzero  elements. 
TABLE  9-3  Maximum  Number  of  Elements  in  Explicit  Transformation  Kernel 

XS     XS(CUB)     XS(PN)     XS(GN) 


GTE 

13 

0 

0 

0 

INVEST 

9 

1 

1 

TANKER 

1 

1 

0 

HFDF 

3 

1 

1 

GAS  PN  A 

1,550 

891 

30 

54 

GAS  PN  C 

18 

5 

0 

0 

GAS  PN  E 

263 

110 

7 

5 

GAS  GN  A 

1,470 

758 

340 

59 

GAS  GN  C 

19 

7 

1 

0 

KELLOGG  2 

6 

2 

0 

0 

KELLOGG  3 

9 

4 

0 

0 

KELLOGG  4 

79 

41 

10 

11 

KELLOGG  5 

1,299 

1,015 

138 

128 

ODS  1 

9 

1 

1 

5 

ODS  3 

0 

0 

0 

0 

TAM  5 

28 

22 

15 

18 

TAM  8 

20 

15 

8 

10 

TAM  12 

38 

42 

16 

21 

PHOENIX  10 

32 

21 

21 

1 

PHOENIX  30 

169 

157 

185 

1 

EA6B 

1,155 

270 

0 

0 

DEC 

218 

39 

49 

46 

AMMO  4H 

235 

92 

0 

0 

BUSCH  4 

10 

5 

0 

0 

BAR 

290 

163 

70 

41 

The  number  of  nonzero  elements  (in  nearest  thousands)  in  the  explicit  transforma- 
tion kernel  A^  gives  some  indication  of  how  much  information  is  not  captured  by 
factorization,  as  well  as  an  idea  of  relative  storage  requirements. 


We  see  that  the  maximum  size  of  the  explicit  transformation  kernel  tends  to  decrease  as  the 
generality  of  the  factorization  increases.  Recalling  the  definition  of  the  explicit  transformation 
kernel. 


[-1 


—    {En  -  E\2F22  F2\) 


this  trend  is  as  we  would  expect.  Each  potentially  binding  explicit  row  which  can  be  converted 
to  a  factored  row  reduces  the  likely  size  of  An].  Also,  the  density  of  the  term  -E\2F22F2\ 
generally  increases  with  the  density  of  F^2  ■  For  F22  fc-by-fc,  the  number  of  nonzeros  in  F2^1 
for  the  GUB  factorization  is  fc,  for  pure  networks  perhaps  as  large  ^-,  and  for  generalized 
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networks  as  large  as  k2.  There  are  some  exceptions  to  this  trend  in  Table  9-3,  especially  in 
model  instances  in  which  the  size  of  the  (GN)  factored  row  set  is  not  significantly  larger  than 
that  of  (PN).  This  is  because  for  a  given  factored  kernel  F22,  the  exact  (PN)  representation  of 
F22    is  generally  more  sparse  than  that  of  the  less  exact  floating-point  representation  of  (CN). 
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It  is  usually  the  case  that  many  constraints  are  not  binding  at  optimality,  as  can  be  seen  in 
Table  9-4. 

TABLE  9-4  Binding  Explicit  Constraints  at  Optimality 

Binding/%      GUB/%         PN/%         GN/% 


GTE 

554/58 

20/02 

20/02 

18/02 

INVEST 

763/57 

201/15 

195/15 

162/12 

TANKER 

51/61 

60/72 

39/47 

9/11 

HFDF 

58/95 

29/48 

29/48 

28/46 

GAS  PN  A 

3,068/45 

1,953/29 

299/04 

349/05 

GAS  PN  C 

2,324/61 

850/22 

68/02 

90/02 

GAS  PN  E 

901/76 

573/48 

90/08 

79/07 

GAS  GN  A 

3,064/45 

1,854/27 

1,155/17 

359/05 

GAS  GN  C 

2,323/61 

861/23 

402/11 

89/02 

KELLOGG  2 

1,950/51 

1,015/27 

92/02 

118/03 

KELLOGG  3 

2,942/51 

1,368/24 

148/03 

183/03 

KELLOGG  4 

4,136/54 

2,491/33 

391/05 

452/06 

KELLOGG  5 

5,270/55 

2,305/24 

592/06 

617/06 

ODS  1 

361/12 

83/03 

76/03 

117/04 

ODS  3 

410/69 

0/00 

0/00 

0/00 

TAM  5 

264/60 

227/52 

168/38 

186/42 

TAM  8 

279/66 

240/57 

142/34 

175/42 

TAM  12 

412/66 

352/53 

205/33 

251/40 

PHOENIX  10 

1,098/68 

1,096/68 

1,090/67 

80/05 

PHOENIX  30 

3,298/77 

3,236/75 

3,239/75 

110/03 

EA6B 

2,939/99 

1,334/45 

26/01 

27/01 

DEC 

1,328/61 

736/34 

726/33 

421/19 

AMMO  4H 

6,686/46 

3,153/22 

13/00 

14/00 

BUSCH  4 

840/67 

481/39 

5/00 

8/01 

BAR 

4,687/63 

3,250/44 

1,895/25 

1,450/19 

Not  all  constraints  are  binding  at  optimality.  The  first  column  lists  the  number 
of  binding  explicit  constraints  and  expresses  this  as  a  percentage  of  all  constraints; 
the  following  columns  display  the  number  of  binding  explicit  constraints  and  their 
corresponding  percentages  under  the  alternate  factorizations. 


A  distinguishing  feature  of  dynamic  factorization  is  the  ability  to  limit  attention  to  binding 
constraints,  handling  binding  factored  constraints  with  great  efficiency,  and  working  with  a 
relatively  small  number  of  binding  explicit  constraints.  Explicit  binding  constraints  on  the 
order  of  a  few  thousand,  or  less,  are  quite  manageable.  This  is  well  beyond  the  size  of  previously 
reported  implementations. 
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10      Conclusions 

Previous  research  by  others  generally  suggests  that  specialized  algorithms  such  as  those  pre- 
sented here  are  useful  only  when  the  factored  structure  completely  dominates.  There  are  even 
reports  of  algorithms  for  solving  problems  having  a  single  unfactored  (explicit)  constraint  (Hultz 
and  Klingman  [1978],  Klingman  and  Russell  [1978]).  When  implementations  have  been  re- 
ported, problem  suites  have  been  limited  to  instances  having  a  very  small  number  of  explicit 
constraints,  typically  in  the  range  from  one  to  twenty  (Chen  and  Saigal  [1977],  Glover,  Karney, 
Klingman  and  Russell  [1978],  Glover  and  Klingman  [1981]).  The  consensus  seems  to  be  that 
such  algorithms  are  quite  delicate,  and  deserve  to  be  viewed  as  specialized  algorithms,  useful 
only  for  solving  very  special  problem  instances. 

We  refute  this  view.  Dynamic  factorization  is  competitive  with  commercial-quality  opti- 
mization systems  on  every  model  instance  we  have  tested. 

The  development  here  stresses  the  similarities  among  the  algorithms  and  the  natural  exten- 
sions leading  from  one  to  the  next.  This  is  in  contrast  to  the  development  reported  for  similar, 
non-dynamic  algorithms  (e.g.,  Dantzig  and  Van  Slyke  [1967],  Klingman  and  Russell  [1978], 
and  Hultz  and  Klingman  [1978])  in  which  the  specifics  of  the  individual  algorithm  obscure  the 
generality  of  the  approach.  The  conceptual  difference  between  our  algorithms  is  seen  to  be 
largely  isolated  to  the  structure  of  a  single  algebraic  entity,  the  factored  kernel.  By  abstracting 
the  structure  of  the  factored  kernel  and  concentrating  on  the  general  algorithm  design,  the 
versatility  and  flexibility  of  this  approach  is  clarified. 

The  algorithmic  development  leads  directly  to  an  implementation.  The  resulting  software 
suite  exhibits  a  "single  system  image".  The  modularity  of  the  algorithm  allows  the  definition 
of  an  "abstract  data  type"  (see,  e.g.,  Aho,  Hopcroft  and  Ullman  [1974])  which  isolates  the  data 
structures  and  update  procedures  for  the  factored  kernel  from  the  rest  of  the  implementation. 
Each  factorization  is  seamlessly  integrated  within  the  system  design. 

The  early  1980s  produced  a  great  deal  of  research  on  automatic  identification  of  special 
structure  in  LP  models  (see,  e.g.,  Gunawardane  and  Schrage  [1977],  Glover  [1980],  Schrage 
[1981],  Brown,  McBride  and  Wood  [1985],  and  Bixby  and  Fourer  [1986]).  We  have  incorporated 
the  most  useful  of  these  ideas  into  our  implementation,  and  we  have  what  we  believe  to  be  the 
first  complete  implementation  which  supports  automatic  identification  of  factored  row  sets.  This 
capability  may  be  used  to  identify  new  factored  structure  or  to  validate  or  augment  a  modeler- 
provided  recommendation.  When  faced  with  the  choice  of  either  solving  an  unfactored  model 
instance  or  automatically  identifying  a  factored  structure  and  then  using  the  corresponding 
solver,  our  results  show  that  the  latter  is  nearly  always  to  be  preferred.  Modelers  have  conducted 
extensive  additional  computational  experimentation  with  the  X-System  not  reported  in  this 
paper.  These  results  suggest  that  in  addition  to  the  quantity  of  factored  rows,  the  quality  of 
these  rows  influences  the  performance  of  factorization  algorithms.  While  not  well  understood, 
it  is  clear  that  the  myopic  approach  of  our  heuristics  is  no  substitute  for  the  modeler's  guidance 
in  identifying  factored  structure. 

Processing  networks  (Koene  [1982])  are  network  models  which  allow  proportional  flow  re- 
strictions on  the  arcs  entering  or  leaving  some  nodes.  One  formulation  of  such  a  model  results 
in  a  pure  or  generalized  network  structure  with  a  set  of  complicating  columns.  Chen  and  En- 
gquist  [1986]  propose  a  primal  partitioning  algorithm  for  solving  processing  network  problems. 
An  alternate  formulation  yields  a  pure  or  generalized  network  structure  with  complicating  rows: 
this  is  precisely  the  structure  dealt  with  here. 

The  multicommodity  capacitated  transshipment  problem  (MCTP)  has  been  the  subject  of 
much  research  over  the  years,  and  a  number  of  specialized  algorithms  have  been  proposed  to 
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solve  it  (see,  e.g.,  Assad  [1978]  or  Kennington  [1978]).  The  usual  MCTP  formulation  is  a  pure 
network  which  each  commodity  uses  independently  in  its  own  flow  model,  but  with  side  con- 
straints on  the  total  common  flow  of  all  commodities  over  some  of  the  network  arcs.  The  side 
constraints  form  a  CUB  row  set,  while  the  rest  of  MCTP  forms  a  pure  network;  either  view 
might  be  preferred  depending  upon  size  of  the  common  network,  the  number  of  side  constraints, 
and  the  number  of  commodities.  In  our  experience,  the  network  factorization  usually  domi- 
nates the  GUB  factorization,  and  the  pure  network  factorization  presented  here  is  a  powerful 
technique  for  solving  MCTPs.  As  an  experiment,  we  customized  our  (PN)  implementation  for 
MCTP  to  exploit  the  special  structure  of  the  explicit  side  constraints.  This  highly-specialized 
implementation  performed  no  better  on  AMMO  4H,  and  we  now  believe  that  this  would  be 
true  for  most  MCTPs. 

There  are  problems  which  would  exhibit  a  large  factored  row  structure  if  not  for  a  set  of 
complicating  columns  (e.g.,  see  Brown,  McBride  and  Wood  [1985]).  One  would  expect  the 
structure  of  the  factored  kernel  to  be  dominated  by  that  of  the  predominant  row  structure, 
with  only  occasional  complications  due  to  the  exceptional  columns.  One  might  allow  for  this 
exceptional  structure  in  the  factored  kernel  by  identifying  it  "on-the-fly"  as  the  algorithm 
progresses,  and  preserving  the  sanctity  of  the  core  factorization.  Though  conceptually  simple, 
some  iterations  of  this  algorithm  would  border  on  the  spectacular.  This  approach  may  be 
thought  of  as  a  hybrid  between  the  dynamic  factorization  developed  here  and  dynamic  basis 
triangulation  methods  (see,  e.g.,  Hellerman  and  Rarick  [1971]  and  [1972],  Saunders  [1976]  and 
McBride  [1980]). 

Dynamic  extrinsic  factorization  is  subsumed  by  the  algorithms  presented  in  this  paper  if 
we  activate  functions  in  the  update  analogous  to  the  secondary  exchanges  now  employed.  Es- 
sentially all  that  has  to  be  done  is  ensure  that  successive  factored  components  retain  their 
stipulated  special  structure.  We  speculate  that  this  will  work  best  in  cases  where  model  struc- 
ture is  amenable,  and  quite  likely  will  require  some  model-specific  customization  to  perform 
well  on  difficult  models.  We  have  limited  our  experimentation  to  those  static  extrinsic  cases 
which  we  believe  to  be  most  generally  useful. 
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